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Spin-dynamics simulations have been used to investigate the dynamic behavior of RbMnF3, 
treating it as a classical Heisenberg antiferromagnet on a simple cubic lattice. Time-evolutions of 
spin configurations were determined numerically from coupled equations of motion for individual 
spins using a new algorithm which is based on Suzuki- Trotter decompositions of exponential 
operators. The dynamic structure factor was calculated from the space- and time-displaced spin- 
spin correlation function. The crossover from hydrodynamic to critical behavior of the dispersion 
curve and spin- wave half- width was studied as the temperature was increased towards the critical 
value. The dynamic critical exponent was estimated to be z = (1.43 ± 0.03), which is slightly 
lower than the dynamic scaling prediction, but in good agreement with a recent experimental 
value. Comparisons are made of both the dispersion curve and the lineshapes obtained from our 
simulations with very recent experimental results for RbMnF3 are presented. 
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§1. Introduction 

RbMnF3 has been the subject of numerous experimen- 
tal and theoretical investigations since it is a good phys- 
ical realization of an isotropic three-dimensional Hi 
berg antiferromagnet. Early experimental studiescr 
showed that the Mn 2+ ions, with spin S = 5/2, form a 
simple cubic lattice structure with a nearest-neighbor ex- 
change constant J exp = (0.58 ± 0.06) meV and a second- 
neighbor interaction constant of less than 0.04 meV [both 
defin ed using the exchange constant to be shown in Eq. 
(3.4b the normalization used here differs from that of 
Ref.y by a factor of two] . Magnetic ordering with anti- 
ferromagnetic alignment of spins occurs below the criti- 
cal temperature T c — 83K. The magnetic anisotropy is 
very low, about 6 x 10~ 6 of the exchange field, .and no 
deviation from cubic symmetry was seen at Tc.B'EP 

Both the static properties and the dynamic response 
of RbMnF3 have been examined through neutron scat- 
tering experiments. Windsor and co-workerst§ looked at 
spin- waves at low temperatures and mapped out_the dis- 
persion curve. The early work of Tucciarone et alu found 
that in the critical region the neutron scattering function 
has a central peak (peak at zero frequency transfer.) and 
a spin-wave peak. Later experiments by Cox et afa ob- 
served a small central peak below T c as well. The more 
recent study by Coldea et aM> also found central peaks 
for T < T c , in agreement with previous work. From 
thp_theoretical side, renormalization-group (RNG) below 
TjilP predicts spin- wave peaks, and a central peak in the 
longitudinal component of the neutron scattering func- 
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tion; however, _at T c both renormalization-group0' and 
mode-couplingO theories predict only the presence of 
a spin-wave peak. The experimentally observed central 
peak is thought to be caused by spin diffusion n rcsulting 
fronuionlinearities in the dynamical equations. 0^ Coldea 
et am 1 also obtained the most precise experimental esti- 
mate of the dynamic critical exponent, z = (1.43±fli)4). 
This is slightly smaller than the predicted valueH of 
z = 1.5 for an isotropic Heisenberg antiferromagnet in 
d — 3 dimensions. 

Extensive Monte Carlo studies measured the static 
properties of classical Heisenberg magnets but could not 
examine the true dynamics of the systems. Several large- 
scale spin-dynamics simulations-have probed the behav- 
ior of various classical systemsjoES' however, there are 
no direct comparisons to experimental data for physi- 
cal systems. In the present work we report large-scale 
simulations of the dynamic behavior of the Heisenberg 
antiferromagnet on a simple cubic lattice, and make di- 
rect comparison with experimental data. 

§2. Model and Methods 

The classical Heisenberg antiferromagnet is defined by 
the Hamiltonian 



n 



j ^ s r 

<rr'> 



(2.1) 



where Sr = {Sr x , Sr v , Sr 2 ) is a three-dimensional classi- 
cal spin of unit length at site r and J > is the antifer- 
romagnetic coupling constant between nearest-neighbor 
pairs of spins. All simulations were performed using 
L x L x L simple cubic lattices with periodic bound- 
ary conditions. The dynamics of the spins are governed 
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by the coupled equations of motion^ 



SrxJ Sr '' 
<rr'> 



(2.2) 



and the time dependence of each spin can be determined 
from the integration of these equations. 

The dynamic structure factor 5(q, to) for momentum 
transfer q and frequency transfer u> can be measured by 
inelastic neutron scattering experiments and is given by 



., r+oo 

S* fc (q, w) = e / e tut C k (R, t) 
R J -°° 



dt 



2tt 



(2.3) 



where R = r — r', C fc (R, t) is the space-displaced, time- 
displaced correlation function, with k — x,y, or z, and 

C fc (R,i) = (S r k (t)Sr' k (0)) - (S r k (t))(S r ' k (0)). (2.4) 

In the case of antiferromagnets, the wave- vectors are 
measured with respect to the (ir, it, it) point which cor- 
responds to the Brillouin zone center. Note that in the 
[1,1,1] and [1,0,0] directions the respective first Brillouin 
zone boundary wave-vectors are (±7r/2, ±7r/2, ±7r/2) 
and (±tt,0,0). 

Using Monte Carlo and spin-dynamics methods JiillL 
we simulated the simple-cubic classical Heisenberg anti- 
ferromagnet with 12 < L < 60 at the critical tempera- 
tureEP T c = 1.442929(77) J as well as below T c . (We use 
units in which Boltzmann's constant ks = 1.) Equilib- 
rium configuratipas were generated using a hybrid Monte 
Carlo methodli-ilLa) and the coupled equations of motion 
were then integrated numerically, using these states as 
initial spin configurations. Numerical integrations were 
performed to a maximum time t max < 1000J -1 , using 
a time step of At. The space-displaced, time-displaced 
correlation function C k (R, t) was computed for time- 
displacements ranging from to t cuto ff and extracted 
from an average over 40 to 80 different time starting 
points, evenly spaced by lOAi. As many as 7000 ini- 
tial configurations were used, although for large lattices 
this was reduced to as few as 400. For L — 24 at 
T = 0.9T C the integration was carried out with a time 
step At = QJ31J -1 using a 4th-order predictor-corrector 
method.oES* For other lattice sizes and temperatures, 
we used a new algorithmic EiS based on 4th-order Suzuki- 
Trotter decompositions of exponential operators, with a 
time step At = 0.2 J -1 . The larger integration time step 
allowed us to extend the maximum integration time to 
much larger values than was previously possible. 

In order to reduce the computer resoupcps-nceded we 
calculated partial spin sums "on the fly")l3>E3 however, 
data could then only be kept for the (g, 0,0), (q,q, 0) 
and (g, q, q) directions with q determined by the periodic 
boundary conditions, 

<?=— , n = ±l,±2,...,±(L-l),L. (2.5) 

Since all three Cartesian spatial directions are equivalent 
by symmetry, results for (q, 0,0), (0, g, 0) and (0,0, q) 
were averaged. Similarly, the same operations carried 
out for the (q, q, 0) and (q, q, q) directions were also av- 
eraged over the equivalent reciprocal lattice directions. 



For the Heisenberg ferromagnet the order parameter 
is the total magnetization, which is a conserved quantity. 
Thus the dynamic structure factor 5(q, oj) can be sepa- 
rated into a component along the axis of the total magne- 
tization (longitudinal component) and a transverse com- 
ponent; however, for the isotropic antiferromagnet con- 
sidered here the order parameter is not conserved and 
the longitudinal and transverse components of S"(q, ui) 
cannot be separated in the simulation. Henceforth we 
will use the term dynamic structure factor to refer to 
the average. 

Two practical limitations on spin-dynamics techniques 
are the finite lattice size and the finite evolution time. 
The finite time cutoff can introduce oscillations in 
S(q, uj), which can be smoothed out by convoluting the 
spin-spin correlation function with a resolution function 
in frequency. In neutron scattering experiments the di- 
vergence of the neutron beam gives rise to an intrin- 
sic Gaussian resolution function in q and u and the 
smoothed dynamic structure factor is 



~^ J— tcutoff Z7T 



R 



. / ^ fc (q,c')e ~^du', (2.6) 
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where <5 W is a parameter characterizing the Gaussian res- 
olution function and has to be chosen properly so that 
effects due to the finite time cutoff can be neglected. The 
momentum dependent susceptibility, xfCl); ^ s given by 



r^ fc (q,^=x 5 fe (q) 
J~oo *™ 



(2.7) 



Finite-size scaling theory0EJ can be used to extract 
the dynamic critical exponent z: the divergence of the 
correlation length £ is limited by L and the dynamic 
finite-size relations are given by 

=G{uL*,qL,5 u L*) (2.* 



xi(q) 



and 



L- z n{qL,S^L 2 



where Lo m is a characteristic frequency, defined as 
Sl (q,w)— = -x£(q). 

Z7T Z 



(2.9) 



(2.10) 



For t cu toff > 400 J^ 1 the oscillations in the dynamic 
structure factor were not very significant. Thus, we first 
estimate the dynamic critical exponent z without using 
a resolution function, i.e. we take d u = 0. In this case 
z can be obtained from the slope of a graph of log(w m ) 
vs log(L) (where Lu m is the characteristic frequency for 
S u = 0) if qL is fixed and L is large enough to be in the 
asymptotic-size regime. 

The effects of the small oscillations in S(q,ui) on the 
dynamic exponent z can be evaluated by repeating the 
analysis using a resolution function so that the function 
Vl^qL.S^L 2 ) in Eq. (|2.9|) is constant if qL and 5ujL z are 
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fixed and 



(2.11) 



Because 8 U depends on z, this exponent had to be de- 
termined iteratively. An initial value was used to 
determine 5 U for different L, Sl (q, to) and ui m were 
computed for different values of L and q with qL held 
fixed (i.e. n is const ant) and a new estimate, was 
extracted from Eq. ( 2.11 ). This procedure is repeated 
until the estimates converge. 

§3. Results 

3.1 Numerical data for S((\,oj) 

For T < T c our results for the dynamic structure factor 
show a spin-wave and a central peak. At low tempera- 
tures the central peak is barely visible and very narrow 
spin wave peaks are the dominant feature, see Fig. 1 
(finite-size effects are evident for n = 1). In Fig. 2 we 
show lineshapes for lattice size L = 60 at T — 0.9T C and 
several q-values in the [100] direction. As q increases, 
the central peak broadens and its relative amplitude in- 
creases. As the temperature is increased, the central 
peak grows, and at T c the central peak is even stronger 
than the spin wave peak. Fig. 3 shows lineshapes for 
L = 60 and q = tt/15 and 37r/10 in the [100] direction. 
Clearly oscillations due to the finite t cuto ff are negligi- 
ble; therefore, in our lineshape analysis we did not use a 
resolution function. For direct comparison with experi- 
mental data we convoluted our results with a Gaussian 



3 
CO 





1 1 

:; n=i 










n=2 








n 


=3 




J | 


t 


n=4 


I 

. — ,=ZJ 




I J 


l. A A 



1 2 3 4 5 

co/J 

Fig. 1. Dynamic structure factor S(Q,w) from our simulations 
for L = 20 at T = 0AT c , q in the [100] direction. The symbols 
represent spin dynamics data and the solid line is a fit with the 
Lorentzian function given in Eq. ( |3.l| ). 



resolution function with the same width as in the exper- 
iment. The structure in the lineshapes discussed here is 
much larger than oul. resolution in frequency Below T c , 
previous theoreticaltliP and experimentalET studies pro- 
vided the comparison of the position and the half-width 
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Dynamic structure factor S(q, u>) for L = 60 at T = 
and q in the [100] direction. The symbols represent spin 
dynamics data and th e so lid line is a fit with the Lorentzian 
function given in Eq. (3.1). 



Fig. 2. 
0.9T t 



of the spin-wave and central peaks by fitting the line- 
shape to a Lorentzian form 



S(q,w) 



ATf 



Us) 2 

(3.1) 

where the first term corresponds to the central peak and 
the last two terms are from the spin-wave creation and 
annihilation peaks at w = ±u> s . 

For T — 0.9T C Lorentzian lineshapes fit our results 
well for small values of q, except for the smallest value, 
q = 2ir/L, in the [100] direction for which finite size ef- 
fects are apparent. For large values of q the Lorentzian 
form given in Eq. ( |3 . l[ ) does not fit the data, especially 
at high frequency. In general, the fitted parameters var- 
ied when different frequency ranges were used in the fit 
by an amount which was often larger than the statistical 
error in the fitted parameters obtained from the fit using 
a single frequency range. Therefore, for T = 0.9T C we es- 
timated the error in the fitted parameters by fitting the 
lineshapes using three different ranges of frequency and 
taking thepaverage. At T c , renormalization-group the- 
ory (RNG)lilP predicts a non-Lorentzian functional form 
for the spin-wave lineshape, which has been used along 
with a Lorentzian central peak to analyze experimental 
data.& ) Since it is more complicated to perform fits to 
this RNG functional form and since the spin- wave peaks 
obtained from the simulations are more pronounced than 
in the experiment, and thus less dependent on the fit- 
ted functional, we have fitted the lineshapes at T c to 



Lorentzians, as given in Eq. (3.1). Although obtaining 
a good fit to our data at T c was more difficult than be- 
low T c , the resulting fits are still reasonable. Unlike for 
T = 0.9T C , at T c the lineshape parameters used in the 
analysis below are the values obtained from the fit to 
only one frequency range, which was the one that gave 
the best fit. The actual error in the fitted parameters at 
T c should be larger (by up to a factor of 5) than the error 
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Fig. 3. Dynamic structure factor S(q, u>) at T = T c : (a) q = n/15 
and (b) 37r/10 in the [100] direction. The symbols represent spin 
dynamics data for L = 60 a nd t he solid line is a fit to the sum 
of Lorentzians given in Eq. ( h.l| ). Error bars arc only shown for 
a few typical points; at high frequency they are of the size of the 
fluctuations in the data. 



q values of the dispersion curve to Eq. (3J!) yielded 
x = 1.017 ± 0.003. If larger values of q are included in 
the fit, the exponent decreased slightly. The sensitivity 
of the fitted exponent to the particular form of the fitted 
function was examined by including a quadratic term, 
i.e., 



A s q x + B s q 2 



(3.3) 



yielding a similar value x = 1.020 ± 0.003. When larger 
values of q were included in the fits, Eq. ( |3.3[ ) tended 
to yield smaller x 2 's per degree of freedom than Eq. 
(3.2). The dispersion curve for T — T c and L — 60 fit- 
ted to Eq. (3.2) yielded an exponent of x = 1.38 ± 0.01 
when the smallest 12 values of q were included in the 
fit. As the larger q were excluded from the fit, the 
exponent increased slightly, tending towards x ~ 1.40. 
When only the smallest few values of q were included in 
the fit, the exponent decreased again, reflecting the fact 
that we probed correlations between spins separated by 
larger distances, i.e. smaller q. This reflected the finite 
size of the lattice (and thus of the correlation length), 
showing that the system is not at criticality. Hence, the 
exponent x decreases towards unity. In contrast, large 
values of q correspond to short distance (in the direct 
lattice space) spin-spin correlations, and the correlation 
length is much larger than the distance probed. Our re- 
sults at T c agree with recent experimenttx which found 
x = 1.43 ± 0.04 when the dispersion curve at T c was 
fitted to a power-law relation of the form given in Eq. 
(U). The solid lines in Fig. 4 are fits to Eq. (|J); in 
general, these fits gave lower values of x 2 P er degree of 
freedom than fits to Eq. (jjT^).,— In the critical region, 
dynamic scaling theory predictsElP that the half-width 
of spin- wave peaks behaves as T2 ~ q 1 ' 5 , whereas for 
the hydrodynamic regime the prediction from hydrody- 
namic theoryEP is r 2 ~ q 2 . The half- width of the spm- 



bars show n in the figures below. Illustrations of the fits 
using Eq. (3.1) to the simulated lineshapes at T = 0.9T C 
and at T c for several values of q are shown in Figs. 2 and 
3. 

We also found very weak peaks in the high frequency 
tail of the spin- wave peaks. Using the spin- wave frequen- 
cies in the [100], [110] and [111] directions we could check 
that the position of these extra peaks corresponded to 
frequencies of two spin- wave addition peaks. These ex- 
tra structures in the lineshapes were particularly visible 
for the smallest q- values. 

Fig. 4 shows how the dispersion curve varies as the 
temperature increases from T — to T c . Although the 
Lorentzian in Eq. (3.1) did not yield good fits to the 



lineshapes for larger values of q, the spin-wave peak po- 
sitions could still be determined relatively accurately and 
the dispersion curve could be measured up to q = ir/2. 
Well below T c , the dispersion relation is linear for small 
q, but as the temperature increases towards T c , the dis- 
persion relation changes gradually to power-law 

uj s = A s q x . (3.2) 

For T = 0.9T C and L = 60 a fit to the smallest five 




q (7i) 

Fig. 4. Spin- wave dispersion relations for T < T c , in the [100] 
direction. The symbols represent spin-wave positions extracted 
from Lorentzian fits to the lineshapes from the simulations, and 
the solid curves are fits o f the dispersion relations at different 
temperatures to Eq. (3.3). 



wave peaks at T = 0.9T C and L = 60 from our Simula- 
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tions is shown in Fig. 5. We observed a crossover from 
T 2 = (0.401 ± 0.004)g 1 - 46±0 06 for larger values of q to 
the behavior T 2 = (0.48±0.02)g 1 - 86±005 for small values 
of q. The behavior for relatively large q agrees with dy- 
namic scaling theory and with recent experiment .cf The 
exponent we obtained by fitting only small values of q 
is close to the hydrodynamic prediction. Thus, the spin- 
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Fig. 5. Log-log graph of the half-width of the spin-wave peak 
extracted from Lorentzian fits to the lineshapes obtained from 
simulations for L = 60 and T = 0.9T C in the [100] direction as a 
function of q. 



wave half-width reflected crossover between two different 
regimes, the critical and the hydrodynamic regions. This 
crossover is similar to the one observed in the dispersion 
curve at T c , discussed above. For T = T c and L = 60 
the spin-wave half-width also had a power-law behav- 
ior which varied from approximately q 12 when the 12 
smallest values of q were included to ~ q 1A when only 
the smallest five wave-vectors were considered. In their 
recent experiment, Coldea et a© found T 2 = Dq 1A1±0 - 05 
for 0.77T C < T < T c , and the coefficient D increased with 
increasing temperature. 

As in the experiments, the dynamic structure factors 
from our simulations had central peaks (u> — 0) for 
T <T C . In contrast, RNG theory predicts a central peak 
in the longitudinal component of the dynamic structure 
factor only below T,J3_and none of the theories predict 
a central peak at T c MB> For T = 0.9T C and L = 60 fits 
of the central peak half-width to a simple power law were 
poor, but a much improved fit was obtained by using the 
function Ti = Ai + B\q Cl , which allows for a non-zero 
central peak width when q vanishes. In these fits the data 
for the smallest possible q, i.e. n — qL /2tt = 1, were ex- 
cluded because of large finite-size effects. The fit includ- 
ing data for q corresponding to n — 2 until n — 7 yielded 
Ai ~ 0.013±0.001, B x ~ 0.120±0.005 and d ~ 2.4±0.2. 
As we systematically included larger values of q in the 
fits, these parameters decreased slightly. At T c we also 
fitte d the central peaks to Lorentzians, according to Eq. 
(3.1); however, these tended to yield curves with smaller 
amplitudes than the data. Since there is no theoreti- 



cal prediction for the central peak, we have also tried to 
fit them with a Gaussian form but the result was much 
worse than with Lorentzians. 

The lattice sizes that we used were all multiples of 
12 so there were certain g-values which were common to 
all L. Lineshapes and spin- wave peak positions could 
be compared for different L at a fixed value of q. At 
T = 0.9T C we saw no significant finite-size effect for L > 
24; however, when we superimposed lineshapes at T c for 
fixed q, and different values of L, finite-size effects were 
noticeable for L — 24. For larger L the lineshapes were 
the same within the error bars. 

The dynamic critical exponent z was extracted from 
finite-size scaling of u) m . From an analysis without res- 
olution function, or equivalently S u — 0, and n = 1,2, 
wc estimated L = 30 to be approximately the onset of 
the asymptotic-size regime and z — 1.45 ±0.01 for n = 1 
and z = 1.42 ±0.01 for n — 2. We also estimated the 
value of z using a resolution function. Several initial val- 
ues of z(°) were used, and in all cases the exponent z 
converged rapidly. Our final estimate for the dynamic 
critical exponent is z = 1.43 ± 0.03. Analyses of the 
characteristic frequency function of L with and 

without a resolution function agreed closely. 

3.2 Comparison with experiment for RbMnF^ 

We now compare our results with the recent neutron 
scattering data of Coldea et alW The Mn 2+ ions in 
RbMnF3 have spin S = 5/2 and interact via a quantum 
Heisenberg Hamiltonian of the form 



H = J e 



<rr'> 



(3.4) 



where Sr are spin operators with magnitude |Sr^| 2 = 
S(S + 1) and the interaction strength between pairs of 
nearest-neighbors was determined experimentallyuP to 
be J exp = (0.58 ± 0.06) meV. In contrast, our simu- 
lations were performed on a classical Hamiltonian; how- 
ever, quantum Heisenberg systems with large spin values 
(S > 2) have been shown to behave as classical Heisen- 
berg systems where the spins are vectors of magnitude 
\/S(S + 1) with the same interaction strength between 
pairs of nearest neighbors as in the quantum system. EJ 
Since our classical spins were vectors of unit length, a 
normalization of the interaction strength J from our sim- 
ulation is needed, i.e. 



J = J exp S{S+l). 



(3.5) 



Although this choice leaves the Hamiltonian unchanged, 
it modifies the equations of motions given above as Eq. 
(2.2). The dynamics of the classical system so defined is 



the same as the quantum system defined by the Hamilto- 
nian in Eq. ( |3.4p if one rescales the time, or equivalently, 
the frequency. We obtain 



u exp = J exp ^S(S+ 1) -, 



(3.6) 



where uj exp is the frequency transfer in the quantum 
system, measured experimentally, and w/J is the fre- 
quency transfer in units of J from our simulations. Note 
that the critical temperature of the classical Heisenberg 



6 



Landau, Tsai, and Bunker 



model, has been determined from Monte Carlo simula- 
tionsEP to be T c = 1.442929(77) J. Using the normal- 
ization for the interaction strength J given in Eq. ( |3.5| ) 
and the experimental value J exp — (6.8±0.6)Kl§ we get 
T c = (85.9 ± 7.6)K which is well within the error bars of 
the experimental value of around 83K. 

Neutron scattering experiments measure the dynamic 
structure factor multiplied by a tempepiture and fre- 
quency dependent population factor jaESES and this fac- 
tor was removed from the experimental data before com- 
paring them with the simulational data. Another com- 
plication in the experiment is the finite divergence of 
the neutron beam which gives rise to a resolution func- 
tion which is usually approximated by a Gaussian in the 
4-dimcnskpal energy and wave-vector space. In the ex- 
periment,™ the measured resolution width along the en- 
ergy axis was 0.25 meV (full-width at half-maximum) for 
incoherent elastic scattering. In order to directly com- 
pare our results with the experiment, we convoluted the 
lineshapes from our simulation with a Gaussian resolu- 
tion function in energy with the experimental value of 
full- width at half-maximum, normalized according to Eq. 
(3.6). The standard deviation 8^ thus obtained for the 
Gaussian resolution function in Eq. (2.6) was 0.0619 in 
units of J. (Convolution of our lineshapes with the ex- 
perimental 3-dimcnsional Gaussian function in the wave- 
vector space had a negligible effect.) 

The experimental dataEP are from constant-^ scans 
with both positive and negative energy transfer. The 
wave- vector transfer Q was measured along the [1,1,1] di- 
rection, around the antiferromagnetic zone center which 
injpur notation is the (-7T, 7r, tt) point. Note that Coldea et 
a& define the wave-vector transfer Q in units such that 
the antiferromagnetic zone center is (0.5, 0.5, 0.5); hence, 
to express their Q in units of A -1 one has to multiply 
it by 2n/a, where a is the cubic lattice parameter ex- 
pressed in A. However, in the simulation direct lattice 
positions are defined in units of the lattice constant a; 
thus we obtain wave-vectors multiplied by the constant 
a. We emphasize that one has to divide the wave- vector 
Q [and also q, see Eq. (f2~5|)] defined in this paper by 
27rJn order to express it in the units used by Coldea et 
al.u In the experiment, measurements were taken for 
wave- vectors Q = (n + q, n + q, n + q), with q = 27r(0.02), 
27r(0.04),..., 27r(0.12), but unfortunately these values of 
q are not accessible in our simulations for the particular 
lattice sizes that we used. Thus, direct comparison of 
the lineshapes from the simulation with the experimen- 
tal ones required interpolation of our results to obtain 
the same q values as the experiment. This was done by 
first fitting our lineshapes with a Lorentzian form, as 
given in Eq. (3.1). Since the parameters B, T2 and u> 8 
obtained from these fits behave as power-laws of q, we 
linearly interpolated the logarithm of these parameters 
as a function of the logarithm of q, to obtain new param- 
eters for the lineshapes corresponding to those values of 
q actually observed in the experiment. We estimated the 
uncertainties from this procedure to be less than five per- 
cent for the parameter J5, less than three percent for the 
spin- wave half- width T2 and the spin-wave position cu s at 
T c , and less than one percent for the spin-wave position 



uj s at T — 0.9T C . Below T c , the parameters A and Ti as- 
sociated with the central peak were linearly interpolated, 
yielding new values with uncertainties of approximately 
five percent. At T c , the parameter A was interpolated in 
the log- log plane (as for B, and oj s discussed above), 
whereas Ti was simply linearly interpolated. The uncer- 
tainties in A and Ti at T c were estimated to be less than 
ten percent. For L = 60, there is one value of q, namely 
q = 27r(0.10), which is accessible to both simulation and 
experiment. This was the only case for which we did not 
have to interpolate in q. 

Our results at the critical temperature can be com- 
pared with the experimental data at the same tempera- 
ture. Below T c , the simulations are mainly for T — 0.9T C 
which unfortunately does not coincide with any temper- 
ature used in the experiment; however, it is very close 
to T = 0.894T C for which experimental results are avail- 
able. To correct for the slight difference, we made a 
linear interpolation in temperature, using our results for 
L = 24 at T = 0.846T C and at T = 0.9T C . We first 
fitted the lineshapes at these two temperatures to the 
form given by Eq. ( |3.l|) , then we linearly interpolated 
the position and the amplitude of the spin-wave peak 
at these temperatures, to obtain the spin-wave position 
and amplitude corresponding to T = 0.894T C . For small 
values of q we found that the frequency and amplitude 
of the spin-wave peak at T = 0.894T C were respectively 
~ 1.5 and 5 percent larger than at T = 0.9T C and this 
difference decreased for larger values of q. 

The intensity of the lineshapes in the neutron scatter- 
ing experiment was measured in counts per 15 seconds. 
For both temperatures T = 0.894T C and T = T c the mea- 
surements for the several wave-vectors were done with 
the same experimental set-up and conditions. Therefore, 
the relative intensities of the lineshapes for the different 
wave- vectors is fixed, and equal for both temperatures. 
The intensity of the lineshapes obtained in the simulation 
had to be normalized to the experimental value; however, 
because the relative intensities for different wave-vectors 
is fixed, we have only one independent normalization fac- 
tor for all the wave- vectors at both temperatures. The 
normalization of the intensity was chosen so that the 
spin-wave peak for T = 0.894T C and q = 27r(0.08) from 
the experiment and the simulation matched. This same 
factor was used to normalize the intensities of the line- 
shapes corresponding to the remaining values of wave- 
vectors at T = 0.894T C , and for all cases at T c . 

The final lineshapes for T = 0.894T C , L = 60, and 
two wave-vectors are shown in Fig. 6 together with ex- 
perimental lineshapes for each case. Fig. 7 shows the 
comparison of the dispersion curve from the simulation 
and the experiment at T = 0.894T C . Good agreement be- 
tween our results and experiment can be seen from either 
the direct comparison of the lineshapes, or the compar- 
ison of the dispersion curve. The lincshape intensities 
from simulation and experiment agree over two orders of 
magnitude, from q = 2tt(0.02) to q = 2tt(0.10). Fig. 8 
shows the comparison of lineshapes from the simulation 
and the experiment for T — T c , L — 60, and two values 
of q. The dispersion curve obtained from the simulations 
at T = T c , shown in Fig. 7, is systematically larger than 
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the experimental values. We emphasize that the error 
bars shown for the dispersion curve obtained from our 
simulations at T c reflect only the statistical errors of a 
best fit of the lineshapes. For each q, this fit was done 
with only one range of frequency; hence errors associ- 
ated with the choice of frequency range and the quality 
of the fit were not taken into account. It is reasonable 
to expect that such sources of error would increase the 
error bars by a factor of 5. From the direct compari- 
son of the simulated and experimental lineshapes at T c 
it is difficult to determine the difference in the spin-wave 
frequencies, because the spin- wave peaks from the ex- 
periment are not very pronounced, and their positions 
have to be extracted from the fits of the lineshapes. As 
mentioned earlier, the experimental data at T c were fit- 
ted to a functional form predicted by RNG theory plus 
a Lorentzian central peak. The quality of the fitting 
can be seen in Fig. 8(b) for q = 27r(0.08), along with 
the RNG component of the fit and the prediction by 
mode-coupling theory. Finally, even though at T c the 
lineshape intensities from the simulations for small fre- 
quency transfer tended to be lower as compared to the 
experiment, the agreement is still reasonably good, con- 
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Fig. 7. Comparison of the dispersion curve at T = T c and T = 
0.894T C obtained from simulations for L = 60 (circle) and the 
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without interpolation to match the q values from the experiment. 
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sidering the variation of the intensities over almost two 
orders of magnitude from q = 27r(0.02) to q = 27r(0.12). 

§4. Conclusion 

We have studied the magnetic excitations and the dy- 
namic critical properties of the classical Heisenberg an- 
tiferromagnet on a simple cubic lattice using large-scale 
spin dynamics simulations. A new 4th-order decomposi- 
tion integration technique as implemented in Ref.^ al- 
lowed us to use a larger time integration step thus ex- 
tend the maximum integration time to much larger val- 
ues than was previously possible. 

Below T Cl the dispersion curves were approximately 
linear for small q. Increasing the temperature towards 
the critical temperature the dispersion curve became a 
power-law, reflecting the crossover from hydrodynamic 
to critical behavior. The spin-wave half-width at T — 
0.9T C also showed a crossover from critical behavior at 
large values of q to hydrodynamic behavior at small val- 
ues of q. The dynamic critical exponent was estimated 
to be z — (1.43 ± 0.03) which is in_agreement with the 
experimental value of Coldea et alW and slightly lower 
than the dynamic scaling prediction. 

We made direct, quantitative comparison of both 
the dispersion curve and the lineshapes obtained from 
our simulations with the recent experimental results by 
Coldea et a& for RbMnF 3 . At T = 0.894T C the agree- 
ment was quite good with the major difference being at 
T c where spin-wave peaks from our simulations tended 
to be at slightly larger frequencies than the experimen- 
tal results. Both at T = 0.894T C and at T c the line- 
shape intensities varied over almost two orders of mag- 
nitude from q — 27r(0.02) to q — 27r(0.10) and there was 
good agreement between the intensities from simulation 
and experiment over the entire range. Thus, the simple 
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isotropic nearest-neighbor classical Heisenberg model de- 
scribes the dynamic behavior of this real magnetic sys- 
tem quite well. 
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